miss_dat <- cope_subset1 %>%

  mutate(across(starts_with("recode_race"), ~ ifelse(. == "Prefer not to answer", NA, .))) %>%

  mutate(b_dem_sex = ifelse(b_dem_sex %in% c("Prefer not to say", "Other"), NA, b_dem_sex)) %>%
  
  mutate(across(starts_with("recode_orient"), ~ ifelse(. == "I do not want to respond", NA, .)))%>% dplyr::select(b_dem_sex,recode_orient,recode_race)

colnames(miss_dat) <- c("Biological sex","Sexual orientation","Race")

vis_miss(miss_dat)

missing_summary <- miss_var_summary(miss_dat)

miss_plot <- ggplot(missing_summary, aes(x = reorder(variable, -n_miss), y = n_miss)) +
  geom_bar(stat = "identity", fill = "lightblue",size=0.1) +
  geom_text(aes(label = paste0(n_miss, " (", round(pct_miss, 2), "%)")), 
            hjust = 1,  vjust = 0.5, color = "black") +
  coord_flip() +
  theme_minimal() +
  labs( #title = "Missing Data by Variable",
       x = "Variable",
       y = "Missing Count (Percentage)") +
  theme(axis.text.y = element_text(size = 10),
        axis.title = element_text(size = 12),
        plot.title = element_text(size = 14, face = "bold")) +
 # coord_cartesian(ylim = c(0, 40))+
  theme_minimal()
miss_plot 

gg_miss_upset(miss_dat)

dat_cc <- cope_subset3 %>% 
   filter(recode_race!="Prefer not to answer") %>%  ## 1468
  filter(b_dem_sex!="Prefer not to say" & b_dem_sex!="Other") %>% ## 1447
  filter(orientation!="I do not want to respond") ## 1441

cope_subset2 <- cope_subset1 %>%

  mutate(across(starts_with("recode_race"), ~ ifelse(. == "Prefer not to answer", NA, .))) %>%

  mutate(b_dem_sex = ifelse(b_dem_sex %in% c("Prefer not to say", "Other"), NA, b_dem_sex)) %>%
  
  mutate(across(starts_with("recode_orient"), ~ ifelse(. == "I do not want to respond", NA, .)))

## use the latent variable genderid3 as a surragate for a series of gener identity related quesitons
cope_subset2$genderid3 <- genderid_class

## use the latent variable cope3 as a surragate for a series of coping strategies related quesitons
cope_subset2$cope3 <- cope_class

cope_subset3 <- cope_subset2 %>%
  mutate(family_num =challenge_num_dat$challenge_num3,
         family_cat = challenge_cat_dat$challenge_cat) %>%
  dplyr::select(c("b_response_id", "condition", "b_cdi_mean", "f1_cdi_mean", "recode_race",
           #"b_dem_gender_agender", "b_dem_gender_not_sure",
           #"b_dem_gender_other_please_specify", "b_dem_gender_androgynous",
           #"b_dem_gender_nonbinary", "b_dem_gender_two_spirited", "b_dem_gender_female_to_male_transgender_ftm",
           #"b_dem_gender_trans_female_trans_feminine", "b_dem_gender_trans_male_trans_masculine",
           #"b_dem_gender_gender_expansive", "b_dem_gender_third_gender",
           #"b_dem_gender_genderqueer", "b_dem_gender_male_to_female_transgender_mtf",
          # "b_dem_gender_man_boy", "b_dem_gender_transgender", "b_dem_gender_woman_girl",
           "b_screener_age", "b_dem_sex", "recode_orient",
            #"b_covid_family_family_did_not_enough_enough_money_for_food",
           #"b_covid_family_family_did_not_have_a_regular_place_to_sleep_or_stay",
           #"b_covid_family_i_could_not_attend_school_in_person", "b_covid_family_i_could_not_attend_school_at_all",
           #"b_covid_family_other", "b_covid_family_family_did_not_have_enough_money_for_gas_transportation",
           #"b_covid_family_family_did_not_have_enough_money_to_pay_rent",
           #"b_covid_family_the_covid_19_pandemic_has_not_affected_me_or_my_family_in_these_ways_in_the_past_2_weeks",
           #"b_covid_cope_1_connecting_with_others", "b_covid_cope_1_including_talking_with_people_you_trust_about_your_concerns_and_how_you_are_feeling",
           #"b_covid_cope_1_contacting_a_healthcare_provider", "b_covid_cope_1_drinking_alcohol",
           #"b_covid_cope_1_smoking_more_cigarettes_or_vaping_more",
           "recode_language",
           "genderid3","family_num","family_cat",
           # "family3",
          "cope3")) %>% 
  mutate(
   # condition = factor(condition),
    recode_race = factor(recode_race),
   # b_screener_age = factor(b_screener_age),
    b_dem_sex = factor(b_dem_sex),
    recode_orient = factor(recode_orient),
    genderid3 = factor(genderid3),
    family_num = factor(family_num),
    family_cat = factor(family_cat),
    cope3 = factor(cope3)
  )

dat_cc <- cope_subset3 %>% 
   filter(!is.na(recode_race)) %>%  ## 1468
  filter(!is.na(b_dem_sex)) %>% ## 1447
  filter(!is.na(recode_orient) )## 1441


dat_cc %>%  datatable(
  rownames = FALSE,
  options = list(
    pageLength = 5))

save(dat_cc,file="data/dat_cc.rds")

table1_dat <- cope_subset3 %>% 

 mutate(genderid3 = dplyr::recode(genderid3,
                                "1" = "Non-binary",
                                "2" = "Women/girls",
                                "3" = "Male/Masculine")
   ) %>% 
     mutate(family_num = dplyr::recode(family_num,
                                "0" = "0",
                                "2" = ">=2",
                                "1" = "1")
   )  %>% 
  mutate(cope3 = dplyr::recode(cope3,
                                "3" = "Positive",
                                "2" = "No action",
                                "1" = "Combined")
   )  %>% 
   filter(!is.na(recode_race)) %>%  ## 1468
  filter(!is.na(b_dem_sex)) %>% ## 1447
  filter(!is.na(recode_orient) )## 1441
  

colnames(table1_dat) <- c("id", "condition", "Baseline CDI mean score(0-2)", "3-month CDI mean score", 
"Race", "Age (yrs)", "Biological sex", "Sexual orientation", 
"Language", "Gender identity", "Number of challenges", "Type of challenges", "Type of coping strategies"
)

basemodel_dat <- dat_abc %>% dplyr::select(-condition,-b_cdi_mean,-b_response_id,-b_screener_age,)
set.seed(1017)

trainIndex <- createDataPartition(basemodel_dat$f1_cdi_mean, p = .8, list = FALSE, times = 1)
trainData <- basemodel_dat[trainIndex, ]
testData <- basemodel_dat[-trainIndex, ]



train_control <- trainControl(method = "LOOCV")

#train_control <- trainControl(method = "cv", number = 3)

# rf_model <- train(f1_cdi_mean ~ ., 
#                   data = trainData, 
#                   method = "rf", 
#                   trControl = train_control)
tune_grid <- expand.grid(alpha = 0,
                         lambda = 10^seq(-3, 1, length = 10))

abc_model <- train(f1_cdi_mean ~ ., 
                  data = trainData, 
                  method = "glmnet", 
                  trControl = train_control,
                  tuneGrid = tune_grid,
                  metric = "RMSE")

valid_pred <- predict(abc_model, newdata = trainData)

test_pred <- predict(abc_model, newdata = testData)

## mse
rmse_f <- function(actual, pred) {
  sqrt(mean((actual - pred) ^ 2))
}


## perf metrics
### rmse
valid_rmse <- rmse_f(trainData$f1_cdi_mean, valid_pred)
test_rmse <- rmse_f(testData$f1_cdi_mean, test_pred)

### calibration slope & at large
valid_cali_m <- lm(trainData$f1_cdi_mean~ valid_pred)
valid_cali_slope <- coef(valid_cali_m)[2]
valid_cali_at_large <- coef(valid_cali_m)[1]

test_cali_m <- lm(testData$f1_cdi_mean~ test_pred)
test_cali_slope <- coef(test_cali_m)[2]
test_cali_at_large <- coef(test_cali_m)[1]

perf_abc <- data.frame(
  Metric = c("RMSE", "Calibration Slope", "Calibration at Large"),
  Validation = c(valid_rmse, valid_cali_slope, valid_cali_at_large),
  Test = c(test_rmse, test_cali_slope, test_cali_at_large)
)

basemodel_dat <- dat_person%>% dplyr::select(-condition,-b_cdi_mean,-b_response_id,-b_screener_age,)
set.seed(1017)

trainIndex <- createDataPartition(basemodel_dat$f1_cdi_mean, p = .8, list = FALSE, times = 1)
trainData <- basemodel_dat[trainIndex, ]
testData <- basemodel_dat[-trainIndex, ]



train_control <- trainControl(method = "LOOCV")

#train_control <- trainControl(method = "cv", number = 3)

# rf_model <- train(f1_cdi_mean ~ ., 
#                   data = trainData, 
#                   method = "rf", 
#                   trControl = train_control)
tune_grid <- expand.grid(alpha = 10^seq(0,1, length = 10),
                         lambda = 10^seq(-3, 1, length = 10))
person_model <- train(f1_cdi_mean ~ ., 
                  data = trainData, 
                  method = "glmnet", 
                  trControl = train_control,
                  tuneGrid = tune_grid,
                  metric = "RMSE")

valid_pred <- predict(person_model, newdata = trainData)

test_pred <- predict(person_model, newdata = testData)

## mse
rmse_f <- function(actual, pred) {
  sqrt(mean((actual - pred) ^ 2))
}


## perf metrics
### rmse
valid_rmse <- rmse_f(trainData$f1_cdi_mean, valid_pred)
test_rmse <- rmse_f(testData$f1_cdi_mean, test_pred)

### calibration slope & at large
valid_cali_m <- lm(trainData$f1_cdi_mean~ valid_pred)
valid_cali_slope <- coef(valid_cali_m)[2]
valid_cali_at_large <- coef(valid_cali_m)[1]

test_cali_m <- lm(testData$f1_cdi_mean~ test_pred)
test_cali_slope <- coef(test_cali_m)[2]
test_cali_at_large <- coef(test_cali_m)[1]

perf_person <- data.frame(
  Metric = c("RMSE", "Calibration Slope", "Calibration at Large"),
  Validation = c(valid_rmse, valid_cali_slope, valid_cali_at_large),
  Test = c(test_rmse, test_cali_slope, test_cali_at_large)
)

abc_model %>% summary()

person_model %>% summary()

rownames(perf_abc) <- c()
rownames(perf_person) <- c()

kable(perf_abc, caption = "Model performance of baseline risk model (project ABC)")

kable(perf_person, caption = "Model performance of baseline risk model (project personality)")

basemodel_dat_abc <- dat_abc %>% dplyr::select(-condition,-b_cdi_mean,-b_response_id,-b_screener_age,)

best_lambda_abc <- abc_model$bestTune$lambda
best_alpha <- abc_model$bestTune$alpha ## fixed 0 (ridge)

covariates_abc <- model.matrix(f1_cdi_mean ~ ., data = basemodel_dat_abc)[, -1]  # Create the model matrix
y_abc <- basemodel_dat_abc$f1_cdi_mean 

abc_base_m <- glmnet::glmnet(covariates_abc, y_abc, alpha = best_alpha, lambda = best_lambda_abc)

coef(abc_base_m)

abc_base_m <- lm(basemodel_dat_abc$f1_cdi_mean~.,data = basemodel_dat_abc)


# sjPlot::tab_model(abc_base_m)

basemodel_dat_person <- dat_person %>% dplyr::select(-condition,-b_cdi_mean,-b_response_id,-b_screener_age,)

best_lambda_person <- person_model$bestTune$lambda
best_alpha <- abc_model$bestTune$alpha ## fixed 0 (ridge)

covariates_person <- model.matrix(f1_cdi_mean ~ ., data = basemodel_dat_person)[, -1]  # Create the model matrix
y_person <- basemodel_dat_person$f1_cdi_mean 

person_base_m <- glmnet::glmnet(covariates_person, y_person, alpha = best_alpha, lambda = best_lambda_person)
coef(person_base_m)

person_base_m <- lm(basemodel_dat_person$f1_cdi_mean~.,data = basemodel_dat_person)

# sjPlot::tab_model(person_base_m)

## LRT for each variables

### define model formulas
full_formula <- f1_cdi_mean ~ recode_race + b_dem_sex + recode_orient + genderid3 + family_cat + cope3

gender_formula <- f1_cdi_mean ~ recode_race + b_dem_sex + recode_orient + family_cat + cope3

race_formula <- f1_cdi_mean ~  b_dem_sex + recode_orient + genderid3 + family_cat + cope3

sex_formula <- f1_cdi_mean ~ recode_race + recode_orient + genderid3 + family_cat + cope3


orient_formula <- f1_cdi_mean ~ recode_race + b_dem_sex + genderid3 + family_cat + cope3

family_cat_formula <- f1_cdi_mean ~ recode_race + b_dem_sex + recode_orient + genderid3  + cope3

cope_formula <- f1_cdi_mean ~ recode_race + b_dem_sex + recode_orient + genderid3 + family_cat



full_person <- lm(as.formula(full_formula),data =basemodel_dat_person )
gender_person <- lm(as.formula(gender_formula),data = basemodel_dat_person)
lrt_result <- lmtest::lrtest(gender_person,full_person)

deviance <- lrt_result$Chisq[2]

# Extract the degrees of freedom
df <- lrt_result$Df[2]

# Extract the p-value
p_value <- lrt_result$`Pr(>Chisq)`[2] %>% round(2)




full_formula <- f1_cdi_mean ~ recode_race + b_dem_sex + recode_orient + genderid3 + family_cat + cope3
gender_formula <- f1_cdi_mean ~ recode_race + b_dem_sex + recode_orient + family_cat + cope3
race_formula <- f1_cdi_mean ~ b_dem_sex + recode_orient + genderid3 + family_cat + cope3
sex_formula <- f1_cdi_mean ~ recode_race + recode_orient + genderid3 + family_cat + cope3
orient_formula <- f1_cdi_mean ~ recode_race + b_dem_sex + genderid3 + family_cat + cope3
family_cat_formula <- f1_cdi_mean ~ recode_race + b_dem_sex + recode_orient + genderid3 + cope3
cope_formula <- f1_cdi_mean ~ recode_race + b_dem_sex + recode_orient + genderid3 + family_cat


full_abc <- lm(as.formula(full_formula), data = basemodel_dat_abc)
full_person <- lm(as.formula(full_formula), data = basemodel_dat_person)

lrt_person_f <- function(full_model, reduced_formula) {
  reduced_model <- lm(as.formula(reduced_formula), data = basemodel_dat_person)
  lrt_result <- lrtest(reduced_model,full_model )

  deviance <- lrt_result$Chisq[2]
  df <- lrt_result$Df[2]
  p_value <- round(lrt_result$`Pr(>Chisq)`[2], 2)

  return(data.frame(Deviance = deviance, DF = df, P_value = p_value))
}

lrt_abc_f <- function(full_model, reduced_formula) {
  reduced_model <- lm(as.formula(reduced_formula), data = basemodel_dat_abc)
  lrt_result <- lrtest(reduced_model,full_model)

  deviance <- lrt_result$Chisq[2]
  df <- lrt_result$Df[2]
  p_value <- round(lrt_result$`Pr(>Chisq)`[2], 2)

  return(data.frame(Deviance = deviance, DF = df, P_value = p_value))
}


reduced_models <- list(
  `Gender identity` = gender_formula,
  `Race` = race_formula,
  `Biological sex` = sex_formula,
  `Sexual Orientation` = orient_formula,
  `Challenges` = family_cat_formula,
  `Coping strategies` = cope_formula
)

## lrt for project ABC
dev_results_abc <- lapply(names(reduced_models), function(var) {
  result <- lrt_abc_f(full_abc, reduced_models[[var]])
  result$Variable <- var
  return(result)
})


Dev_df_abc <- bind_rows(dev_results_abc) %>%
  dplyr::select(Variable,everything())

kable(Dev_df_abc)

## perform lrt (person)
dev_results_person <- lapply(names(reduced_models), function(var) {
  result <- lrt_person_f(full_person, reduced_models[[var]])
  result$Variable <- var
  return(result)
})




Dev_df_person <- bind_rows(dev_results_person) %>%
  dplyr::select(Variable,everything())

kable(Dev_df_person)

abc_main_dat <- dat_abc %>%  dplyr::select(f1_cdi_mean,condition,b_cdi_mean)

abc_main_m <- lm(f1_cdi_mean~ b_cdi_mean + condition,data=abc_main_dat)

person_main_dat <- dat_person %>%  dplyr::select(f1_cdi_mean,condition,b_cdi_mean)

person_main_m <- lm(f1_cdi_mean~ b_cdi_mean + condition,data=person_main_dat)

abc_main_m %>% summary()

tab_model(abc_main_m)

person_main_m %>% summary()

tab_model(person_main_m)

lp_abc <- predict(abc_base_m, dat_abc)
lp_person <- predict(person_base_m, dat_person)

abc_hte_m <- lm(f1_cdi_mean~ b_cdi_mean + condition*lp_abc,data=dat_abc)
person_hte_m <- lm(f1_cdi_mean~ b_cdi_mean + condition*lp_person,data=dat_person)

## compute average treatment effect
#ate_dat <- dat_abc

# ate_dat$condition = 'Placebo Control'
# 
# abc_ate_ctrl <- predict(abc_main_m,ate_dat)
# 
# ate_dat$condition  = 'Project ABC'
# 
# abc_ate_trt <- predict(abc_main_m,ate_dat)
# 
# abc_ate_trt-abc_ate_ctrl

abc_ate <- coef(abc_main_m)[3]
person_ate <- coef(person_main_m)[3]

## cATE/HTE
## project ABC vs control
hte_dat <- dat_abc

hte_dat$condition = 'Placebo Control'
abc_hte_ctrl <- predict(abc_hte_m,hte_dat)

hte_dat$condition  = 'Project ABC'
abc_hte_trt <- predict(abc_hte_m,hte_dat)

abc_hte <- abc_hte_trt-abc_hte_ctrl

## hte personality vs control
hte_dat <- dat_person

hte_dat$condition = 'Placebo Control'
person_hte_ctrl <- predict(person_hte_m,hte_dat)

hte_dat$condition  = 'Project Personality'
person_hte_trt <- predict(person_hte_m,hte_dat)
person_hte <- person_hte_trt-person_hte_ctrl


point_df <- data.frame(Comparsion= c("Project ABC vs Control","Project personality vs Control" ),
                       ATE = c(abc_ate,person_ate),
                       SE = c(summary(abc_main_m)$coefficients[, "Std. Error"][3],
                              summary(person_main_m)$coefficients[, "Std. Error"][3])
                       )

kable(point_df, caption = "summary of avearged treatment effect")

## calibration plot
abc_cal_dat <- data.frame(abc_ate,
                      abc_hte,
                      lp_abc)

abc_cal_dat$quantile_grp <- cut(abc_cal_dat$lp_abc,
                           breaks = quantile(abc_cal_dat$lp_abc, probs = seq(0, 1, by = 0.2)),
                           include.lowest = TRUE,
                           labels = FALSE)

## compute the averaged hte effect in each quantile risk grp
abc_avg_hte <- abc_cal_dat %>%
  group_by(quantile_grp) %>%
  summarise(abc_avg_hte = mean(abc_hte, na.rm = TRUE))

## replciate for project personality vs control

## calibration plot
person_cal_dat <- data.frame(person_ate,
                          person_hte,
                          lp_person)

person_cal_dat$quantile_grp <- cut(person_cal_dat$lp_person,
                               breaks = quantile(person_cal_dat$lp_person, probs = seq(0, 1, by = 0.2)),
                               include.lowest = TRUE,
                               labels = FALSE)

## compute the averaged hte effect in each quantile risk grp
person_avg_hte <- person_cal_dat %>%
  group_by(quantile_grp) %>%
  summarise(person_avg_hte = mean(person_hte, na.rm = TRUE))

### this is not correct, should boot from the very beginning
## compute bootstrapped CI
# set.seed(1017)
# person_boot_res <- NULL
# for (i in 1:1000){
#   boot.idx <- sample(1:dim(person_cal_dat)[1], size = dim(person_cal_dat)[1], replace = T)
#   boot.data <- person_cal_dat[boot.idx,]
#   boot.data$quantile_grp <- cut(boot.data$lp_person,
#                                   breaks = quantile(boot.data$lp_person, probs = seq(0, 1, by = 0.2)),
#                                   include.lowest = TRUE,
#                                   labels = FALSE)
#   
#   person_avg_hte <- boot.data %>%
#     group_by(quantile_grp) %>%
#     summarise(person_avg_hte = mean(person_hte, na.rm = TRUE))
#   
#   person_boot_res <- rbind(person_boot_res,person_avg_hte)
# }
# 
# person_boot_cali_ci <- person_boot_res %>%
#   group_by(quantile_grp) %>%
#   summarise(
#     lower_quantile = round(quantile(person_avg_hte, probs = 0.025), 4),
#     upper_quantile = round(quantile(person_avg_hte, probs = 0.975), 4)
#   )
# 
# 
# person_cali_result <- left_join(person_avg_hte,person_boot_cali_ci)

abc_cal_dat$lp_person %>% summary()

person_cal_dat$lp_person %>% summary()

## define bootstrap function
boot_f <- function(mydat=NULL,seed=1017){
  set.seed(seed)
  abc_boot_res <- NULL
  for (i in 1:1000){
    boot.idx <- sample(1:dim(mydat)[1], size = dim(mydat)[1], replace = T)
    boot.data <- mydat[boot.idx,]
    
    abc_main_dat <- boot.data %>%  dplyr::select(f1_cdi_mean,condition,b_cdi_mean)
    
    abc_main_m <- lm(f1_cdi_mean~ b_cdi_mean + condition,data=abc_main_dat)
    
    lp_abc <- predict(abc_main_m,boot.data)
    
    abc_ate <- coef(abc_main_m)[3]
    
    abc_hte_m <- lm(f1_cdi_mean~ b_cdi_mean + condition*lp_abc,data=boot.data)
    
    ## cATE/HTE
    ## project ABC vs control
    hte_dat <- boot.data
    
    trt_levels <- unique(boot.data$condition) %>% as.character()
    
    trt <- trt_levels[!(trt_levels %in% 'Placebo Control')]
    
    hte_dat$condition = 'Placebo Control'
    abc_hte_ctrl <- predict(abc_hte_m,hte_dat)
    
    hte_dat$condition  = trt 
    abc_hte_trt <- predict(abc_hte_m,hte_dat)
    
    abc_hte <- abc_hte_trt-abc_hte_ctrl
    
    
    
    ## calibration plot
    abc_cal_dat <- data.frame(abc_ate,
                              abc_hte,
                              lp_abc)
    
    abc_cal_dat$quantile_grp <- cut(abc_cal_dat$lp_abc,
                                    breaks = quantile(abc_cal_dat$lp_abc, probs = seq(0, 1, by = 0.2)),
                                    include.lowest = TRUE,
                                    labels = FALSE)
    
    ## compute the averaged hte effect in each quantile risk grp
    avg_hte <- abc_cal_dat %>%
      group_by(quantile_grp) %>%
      summarise(avg_hte = mean(abc_hte, na.rm = TRUE))
    
    abc_boot_res <- rbind(abc_boot_res,avg_hte)
  }
  return(abc_boot_res)
}

## compute bootstrapped CI
abc_boot_res <- boot_f(mydat = dat_abc)
abc_boot_cali_ci <- abc_boot_res %>%
  group_by(quantile_grp) %>%
  summarise(
    lower_quantile = round(quantile(avg_hte, probs = 0.025), 4),
    upper_quantile = round(quantile(avg_hte, probs = 0.975), 4)
  )
abc_cali_result <- left_join(abc_avg_hte,abc_boot_cali_ci)

## for personality

person_boot_res <- boot_f(mydat = dat_person)

person_boot_cali_ci <- person_boot_res %>%
  group_by(quantile_grp) %>%
  summarise(
    lower_quantile = round(quantile(avg_hte, probs = 0.025), 4),
    upper_quantile = round(quantile(avg_hte, probs = 0.975), 4)
  )


person_cali_result <- left_join(person_avg_hte,person_boot_cali_ci)

## identify the bounds for cali plots

max_cali_y <- max(abc_cali_result$upper_quantile,person_cali_result$upper_quantile) 

min_cali_y <- min(abc_cali_result$lower_quantile,person_cali_result$lower_quantile) 

## check if the ylims are correct 
### I used -0.2 and 0.5 range for the y axis in the calibration plots

logic1 <- min_cali_y <=  0.05 & min_cali_y >=-0.2

logic2 <- max_cali_y <=  0.05 & max_cali_y >=-0.2

if(!logic1 & ! logic2){
  print(paste0("please make sure the range of y axis in the following calibration plot is [", min_cali_y, ",",max_cali_y,']'))
}else{logic_3 = logic1+logic2}

abc_cali_plot <- ggplot(abc_cali_result, aes(x = as.factor(quantile_grp), y = abc_avg_hte)) +
  geom_point(size = 2.5) +
  # geom_line(aes(group = 1), color = "blue") +
  geom_errorbar(aes(ymin = lower_quantile, ymax = upper_quantile), width = 0.05, color = "black") +
  labs(x = "",
       y = "Mean heter_mean_diff"
       # title = ""
  ) +
  geom_hline(yintercept = abc_ate, linetype = "dashed", color = "grey")+
  theme_minimal()+
  scale_y_continuous(
    limits = c(-0.2, 0.05), 
    breaks = seq(-0.2, 0.05, by = 0.05)
  )+
 # ylim(c(min_cali_y,max_cali_y))+
  ylab("Mean difference")+
  # ylab(expression(atop("Risk Difference, %", atop("Harm" %<-% phantom("     ") %->% "Benefit")))) +  # Custom y-axis label
  theme(axis.title.y = element_text(angle = 90, vjust = 0.5, size = 14))

abc_lp_plot <- ggplot(abc_cal_dat, aes(x = lp_abc)) +
  geom_histogram(binwidth = 0.02, fill = "black") +
  labs(x = "Predicted CDI mean score using baseline covariates", y = NULL) +
  theme_minimal()

abc_combo_plot <- cowplot::plot_grid(
  abc_cali_plot,
  abc_lp_plot,
  ncol = 1,
  align = "v",
  rel_heights = c(3, 1)  
)



# Display the combined plot
abc_combo_plot

abc_cal_dat$lp_abc %>% summary()

person_cali_plot <- ggplot(person_cali_result, aes(x = as.factor(quantile_grp), y = person_avg_hte)) +
  geom_point(size = 2.5) +
  # geom_line(aes(group = 1), color = "blue") +
  geom_errorbar(aes(ymin = lower_quantile, ymax = upper_quantile), width = 0.05, color = "black") +
  labs(x = "",
       y = "Mean heter_mean_diff"
       # title = ""
       ) +
  geom_hline(yintercept = person_ate, linetype = "dashed", color = "grey")+
  theme_minimal()+
  scale_y_continuous(
    limits = c(-0.2, 0.05), 
    breaks = seq(-0.2, 0.05, by = 0.05)
  )+
  ylab("Mean difference")+
 # ylab(expression(atop("Risk Difference, %", atop("Harm" %<-% phantom("     ") %->% "Benefit")))) +  # Custom y-axis label
  theme(axis.title.y = element_text(angle = 90, vjust = 0.5, size = 14))

person_lp_plot <- ggplot(person_cal_dat, aes(x = lp_person)) +
  geom_histogram(binwidth = 0.01, fill = "black") +
  labs(x = "Predicted CDI mean score using baseline covariates", y = NULL) +
  theme_minimal()

person_combo_plot <- cowplot::plot_grid(
  person_cali_plot,
  person_lp_plot,
  ncol = 1,
  align = "v",
  rel_heights = c(3, 1)  
)

person_combo_plot

person_cal_dat$lp_person %>% summary()

cal_plots <- ggarrange(abc_combo_plot, person_combo_plot, ncol = 2, nrow = 1)
cal_plots

ggsave("figures/cal_plots.png", width = 10, height = 7, dpi = 600)

## combined plot for report


abc_lp_plot <- ggplot(abc_cal_dat, aes(x = lp_abc)) +
  geom_histogram(binwidth = 0.03, fill = "#1f78b4",alpha=0.8) +
  labs(x = "Predicted CDI-SF mean score (Project ABC)", y = NULL) +
  theme_minimal()

person_lp_plot <- ggplot(person_cal_dat, aes(x = lp_person)) +
  geom_histogram(binwidth = 0.03, fill = "#8B0000",alpha=0.8) +
  labs(x = "Predicted CDI-SF mean score (Project Personality)", y = NULL) +
  theme_minimal()


cal_plots_report <- ggarrange(abc_lp_plot, person_lp_plot, ncol = 2, nrow = 1)
cal_plots_report

ggsave("figures/cal_plots_report.png", width = 10, height = 7, dpi = 600)


abc_cal_dat_report <- abc_cal_dat

person_cal_dat_report <- person_cal_dat


abc_cal_dat_report$grp ="Project ABC"

person_cal_dat_report$grp = "Project Personality"

colnames(person_cal_dat_report) = colnames(abc_cal_dat_report)

combined_cal_dat_report<- rbind(abc_cal_dat_report,person_cal_dat_report)


combined_hist_plot <-ggplot(combined_cal_dat_report, aes(x = lp_abc, fill = grp)) +
  geom_histogram(binwidth = 0.02, alpha = 0.8, position = "identity") +
  labs(x = "Predicted 3-month CDI-SF mean score using baseline prediction models", y = "Frequency", fill = "Project") +
  theme_minimal() +
  scale_fill_manual(values = c("Project ABC" = "#1f78b4", "Project Personality" = "#8B0000"))

combined_hist_plot

ggsave("figures/combined_hist_plot.png", width = 10, height = 7, dpi = 600)